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Spectral-Domain Computation of Fields Radiated 
by Sources in Non-Birefringent Anisotropic Media 

Kamalesh Sainath and Fernando L. Teixeira 


Abstract —We derive the key expressions to robustly address 
the eigenfunction expansion-based analysis of electromagnetic 
(EM) fields produced by current sources within planar non- 
birefringent anisotropic medium (NBAM) layers. In NBAM, 
the highly symmetric permeability and permittivity tensors can 
induce directionally-dependent, but polarization independent, 
propagation properties supporting “degenerate” characteristic 
polarizations, i.e. four linearly-independent eigenvectors asso¬ 
ciated with only two (rather than four) unique, non-defective 
eigenvalues. We first explain problems that can arise when the 
source(s) specifically reside within NBAM planar layers when 
using canonical field expressions. To remedy these problems, 
we exhibit alternative spectral-domain field expressions, im¬ 
mune to such problems, that form the foundation for a robust 
eigenfunction expansion-based analysis of time-harmonic EM 
radiation and scattering within such type of planar-layered 
media. Numerical results demonstrate the high accuracy and 
stability achievable using this algorithm. 

Index Terms —Stratified media; transformation optics. 

I. Introduction 

Environments with (locally) planar-layered profiles are en¬ 
countered in diverse applications such as geophysical explo¬ 
ration, ground penetrating radar, conformal antenna design, 
and so on 0 0 To facilitate electromagnetic (EM) radiation 
analysis in such environments, eigenfunction (plane wave) ex¬ 
pansions (PWE) have long been used because of their relative 
computational efficiency versus brute-force numerical methods 
such as finite difference and finite element methods. More¬ 
over, PWE can accommodate linear, but otherwise arbitrary 
anisotropic layers characterized by arbitrary (diagonalizable) 
3x3 material tensors jl]] . This proves useful when rigorously 
modeling planar media simultaneously exhibiting both electri¬ 
cal and magnetic anisotropy, such as (i) isoimpedance beam- 
shifting devices and (to facilitate proximal antenna place¬ 
ment) ground-plane-coating slabs systematically designed via 
transformation optics (T.O.) techniques 0. 0. (ii) more 
practically realizable (albeit not necessarily isoimpedance) 
approximations to T.O.-inspired media such as metamaterial- 
based thin, wide-angle, and polarization-robust absorbers to 
facilitate (for example) radar cross section control |5|, as 
well as (iii) numerous other media such as certain types of 
liquid crystals, elastic media subject to small deformations, 
and superconductors at high temperatures (6). These named, 
amongst other, modeling scenarios share in common the 
potential presence of a particular class of anisotropic media in 
which the magnetic permeability (/x r ) and electric permittivity 
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(e r ) tensor properties are “matched” to each other and hence 
together define media supporting four “degenerate” plane wave 
eigenfunctions that, while possessing four linearly independent 
field polarization states (eigenvectors) as usual, share only two 
unique (albeit, critically still, non-defective) eigenvalues q. 
Alternatively stated, propagation characteristics within such 
media are still (in general) dependent on propagation direction 
but independent of polarization, eliminating “double refrac¬ 
tion” (“birefringence”) effects 0 0 Hence our proposed 
moniker “Non-Birefringent Anisotropic Medium” (NBAM), 
rather than the “pseudo-isotropic” moniker (6). 

From an analytical standpoint, said PWE constitute spectral 
integrals exactly quantifying the radiated fields 0. Except for 
some very simple cases however, these expansions must almost 
always be evaluated by means of numerical quadratures or 
cubatures, whose robust computation (with respect to varying 
source and layer properties) is far from trivial and requires 
careful choice of appropriate quadrature rules, complex-plane 
integration contours, etc. to mitigate discretization and trunca¬ 
tion errors as well as accelerate convergence 0 0 0 In ad¬ 
dition to such considerations of primarily numerical character, 
a distinct problem occurs, due to said eigenvalue degeneracy, 
when sources radiate within NBAM layers. Indeed, this case 
requires proper analytical “pre-treatment” of the fundamental 
spectral-domain field expressions to avoid two sources of 
“breakdown”: (i) Numerically unstable calculations (namely, 
divisions by zero) during the computation chain, as well as 
(ii) Corruption of the correct form of the eigenfunctions, viz. 
zexp[ik z z\ instead of the proper form exp [ik z z\, the former 
resulting from a naive, “blanket” application of Cauchy’s 
integral theorem to the canonical field expressions 0 (ig. 

To this end, we first show the key results detailing the de¬ 
generate “direct” (i.e., homogeneous medium) radiated fields 
in the “principal material basis” (PMB) representation with 
respect to which the material tensors are assumed simultane¬ 
ously diagonalized by an orthogonal basis Q Q Subsequently, 
we transform these PMB expressions to the Cartesian basis 
(the PWE’s employed basis). Finally, we employ a robust, 
numerically-stable NBAM polarization decomposition scheme 
to obtain the Cartesian-basis direct field polarization ampli¬ 
tudes. The two-dimensional (2-D) Fourier integral-based PWE 
algorithm, resulting from implanting these derived field ex¬ 
pressions into an otherwise highly robust PWE algorithm [2], 
comprises this paper’s central contribution. 


^ote: The material tensor eigenvectors {vi,V 2 ,V 3 } are not to be con¬ 
fused with the field polarization eigenvectors. 
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II. Problem Statement 


We assume the exp — iut convention in what follows. 
Within a homogeneous medium of material properties 
{e r ,/x r }, the electric field £( r) radiated by electric (J) and 
(equivalent) magnetic (AT) current sources satisfied 


A(-) = V x fi r 1 • V x (•) - k\e r 

•(•) 

(III) 

A{£): 

ik 0 r] 0 J - V x fi r 1 • AT 

(II.2) 

and can be expressed via a 3-D Fourier integral over the field’s 
plane wave constituents {E(k)e zk r }0 

A 1 = Adj (A) /Det (a) 


(H. 3) 

E(k) = A 1 • 

ikoyoJ - V x fi r 1 • M 


(II-4) 
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where, anticipating planar layering along z, the k z spectral 
integral is “analytically” evaluated for every (k x ,k y ) doublet 
manifest in the (typically numerically) evaluated outer 2-D 
Fourier integral. That is to say, by “analytically” evaluated 
we mean that the general (symbolic) closed-form solution of 
the k z integral for arbitrary ( k x ,k y ) doublet, obtained by 
equivalently viewing the k z real-axis integral as a contour 
integral evaluated using Jordan’s Lemma and residue calculus, 
is well-known and can be numerically evaluated at the (k x , k y ) 
doublets 0 ED In particular, analytical evaluation of the k z 
integral yields the “direct” field £ d (r) |1|: 


+oo r 


£ d (v) = 


(2 tt) 2 


t(Az) Y a d e n e fk ™ Az + 


n =1 


i(-Az) Y 


~ Ak nz Az 


n =3 


Jk x Ax+ik y Ay dkxdk (n6) 


where a^(k x ,k y ) is the (source dependent) direct field ampli¬ 
tude of the nth polarization, while e n (k x , k y ) and k nz {k x , k y ) 
are (resp.) the electric field eigenvector (i.e., polarization state) 
and eigenvalue of the nth mode (n = 1,2, 3,4) [ 11. Modes 
labeled with n = 1,2 correspond to up-going polarizations, 
and similarly for down-going modes (n = 3, 4)|j 


2 ko = Wyjn oeo, eo, UO’ Vo — \Zuo/ e O’ e r , and fX r are the vacuum 
wave number, vacuum permittivity, vacuum permeability, vacuum plane wave 
impedance, NBAM relative permittivity tensor, and NBAM relative perme¬ 
ability tensor, respectively. An infinitesimal point/Hertzian dipole current 
resides at r' = (x',y',z'), the observation point resides at r = (x,y,z), 
Ar = r — r' = (Ax, Ay, Az), u(-) denotes the Heaviside step function, 
and k = (k x ,k y ,k z ) denotes the wave vector. Furthermore, fy = pr 1 
and d 0 = k^e zz (r xy r yx - r xx r yy ), where 7 ts = t • tv • s (7 = r, e; 
t, s = x,y,z). All derivations are performed for the electric field, but duality 
in Maxwell’s Equations makes immediate the magnetic field solution. Finally, 
a tilde over variables denotes they are Fourier/wave-number domain quantities. 

3 Adj(-) and Det(-) denote the adjugate and determinant of said argument, 
respectively. Det(A)= do(k z — k\ z )(k z — k 2 z )(k z — k 3z )(k z — k^z), where 
{knz} are the eigenvalues (i.e., longitudinal [ z ] propagation constants). 

4 Please see 111, |11| for other relevant layered-medium expressions. 


The problem with the canonical numerical implementa¬ 
tion of this residue calculus approach lies in its tacit as¬ 
sumption of non-degeneracy (distinctness) in the eigenval¬ 
ues {kiz, k 2z , k 3z , k^z}, which does not hold for NBAM 
media. As an illustration of the polarization-independent 
dispersion behavior of NBAM, consider the dispersion 
relations of a uniaxial-anisotropic medium slab {e r = 
Diag [a, a, b] , p, r = Diag [c, c, d \} (k 2 3 4 = k 2 x + k 2 ) |7], [11]: 

k\ z = [klac- (c/d)^] 1 ^, k 2z = [k%ac - (a/b)k 2 \ , 

k 3z = -ki z , and k± z = -k 2z . Setting {a = y 2 c, b = y 2 d} 

(:y is an arbitrary, non-zero multiplicative constant) renders 
k+ = k\ z = k 2z and k~ = k 3z = k± z , demonstrating the 
plane wave propagation direction dependent , but polarization 
independent , dispersion characteristics of uniaxial NBAM 
This conclusion applies also for more general uniaxial NBAM 
material tensors possessing PMB rotated with respect to the 
Cartesian basis ( 6 ). Similarly, for biaxial NBAM with PMB- 
expressed tensors {/x^ mb = Diag [a, 6 , c ], e^ mb = y 2 /xP mb }, 
the polarization-independent dispersion relations are: 

KT b = [(ykofab - 0 a/c)kl - (b/c)k 2 y ] 1/2 (0.7) 

klT h = ~ [(ykofab - (a/c)k 2 x - (i b/c)k 2 y } 1/2 (II.8) 

with kf = fcL mb = k\T h and k~ = fcr b = kft- 

Now, the two-fold degenerate eigenvalue kf has asso¬ 
ciated with it two linearly independent field polarizations 
describing up-going plane waves ( 6 ); this holds likewise for 
the two down-going polarizations with common eigenvalue 
k~. Mathematically speaking, the eigenvalues {kf,k~} are 
each twice-repeating (i.e., algebraic multiplicity of two) but 
have associated with each of them two linearly independent 
eigenvectors (i.e., geometric multiplicity of two), making them 
non-defective and rendering the four NBAM polarization states 
suitable as a local EM field basis within NBAM layers ©• 
Despite the existence of four linearly independent eigen¬ 
vectors, it is worthwhile to further exhibit the key results 
of the systematic analytical treatment of the two fictitious 

double-poles of A to render numerical PWE-based EM field 
evaluation robust to the two said sources of “breakdown”; this 
treatment is performed in the next section. 

Let us first make two preliminary remarks, however. First, 
assume that the source-containing layer is a biaxial NBAM 
with pf nih — Diag [a, 6, c) and eP mb = ^ 2 /JP mb . Second, 
the orthogonal matrix U = [vi w 2 V 3 ] transforms vectors 
between the PMB and Cartesian basis. For example, the 
relationship between the nth PMB eigenmode wave vector 
kP mb = (fcP“ b ,fcP” b ,A;P” b ) and the (assumed available]^} 
nth Cartesian-basis wave vector k n = (k x ,k y ,k nz ) writes as 

kr b = U- 1 • k n . 

III. Direct Electric Field Radiated within NBAM 

The (Cartesian basis) Fourier domain representation of the 
electric field, radiated in a homogeneous NBAM, writes as 

5 The Cartesian basis wave vectors and polarization eigenvectors are as¬ 
sumed available (e.g., via the state matrix method (TT|). Indeed, recall that the 
operations discussed herein are performed within the backdrop of numerical 
2-D Fourier integral evaluations 0. 
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E = —A • V x fj, r 1 • M for a (equivalent) magnetic 

current source or E = ikopoA • J for an electric current 
source. These two equations, moreover, hold equally when re¬ 
represented in the NBAM’s PMB (i.e., adding “pmb” super¬ 
script to all quantities), which is what we will employ. Indeed, 


the components {A mw } (ra, w = 1, 2, 3) of A (•) write 
as (A rnw — A wrn , and k = k pmb /fco): 

B = -cy 2 kl (k 2 z - [aby 2 - (< a/c)k 2 x - (6/c)fc^]) (III.l) 
A n = (k 2 x - bey 2 ) /B, A 12 = kJ y /B (ffl.2) 

A 13 = kJ z /B, A 22 = (k 2 - acy 2 ) /B (III.3) 

A 2 3 = k y k z /B, A 33 = (, k 2 - aby 2 ) /B (III.4) 

~ — l,pmb . 

while the components of —A • V pmb x ^“ 1,pmb (-) 

mw } write as ( A mw = -A wm ): 

B' = B/(f), A 12 = —ick z /B' (ID.5) 

A 13 = ibky/B 1 , A 23 = —iak x /B' (III. 6 ) 


The expressions within Eqns. ( |III.l| )-( |nL4l ) describe the elec- 
trie field from an electric current source while the expressions 
within Eqns. ( |III.5| )-( |nL6| ) describe the electric field from an 
(equivalent) magnetic current source. Duality in Maxwell’s 
Equations makes immediate the magnetic field results. 

Next the PMB electric field E pmh (k x , k y \ z, z r ), after re¬ 
expressing Eqns. in terms of {k x ,k y ,k z } to 

identify the k z (rather than & pmb ) eigenvalues {k nz } (using 
the relation k = U-k pmb ) as well as “analytically” performing 
the k z contour integral, can be decomposed into a linear 
combination of the degenerate up-going modes {e^ mb ,e 2 mb | 
(for 2 > z f ) or down-going modes {e^ 11113 , e^™ 13 } (for z < z ') |j 
For an electric source, we have for e ±,pmb : 


=b 2tt i 


ik 0 r]o (k z - kf^j A • J pmb 


(HI.7) 


k z =ki 


and similarly for a (equivalent) magnetic source upon replac- 

ing ikorjo A • J pmb with - A • V pmb x //- 1 ’ pmb • 

M p mb i n Eqn. ( |III.7| ). Next, the degenerate PMB modal 
electric fields are re-expressed in the Cartesian basis (e ± = 
U • e ±,pmb ) from which the Cartesian-basis direct field modal 
amplitudes {af , af, af, af} can be robustly extracted using the 
polarization decomposition method proposed previously for 
sources radiating within isotropic layers [ 11: 


r~ 


-1 

r~+-i 

&xl 

&x2 


e' 

c X 

Cyl 

&y2_ 


e + 

L y J 



&x3 

£X4: 

-1 

V 

_&y3 

e y 4_ 


(ifi; 


where {e xn , e yn , e zn } are the x , y , and z components of 
the (cartesian basis) NBAM’s nth electric field eigenvector 
e n . Moreover, if the above-inverted matrices are suspected 
(with respect to, say, the euclidean matrix norm measure) of 
being ill-conditioned, one can always utilize instead say the y 
and z, or alternatively the x and z, components of the field 


6 When z = z', assuming the source does not lie exactly at a planar material 
interface, one can write the direct fields as a linear combination of either the 
up-going or down-going modes since both combinations lead to identical field 
results (save at r') on the plane z = z ' & 0 - 


eigenvectors (TJ. Indeed, this decomposition procedure is well- 
defined due to the non-defective nature of the eigenvalues, 
and hence linear independence between the four NBAM field 
eigenvectors {e n } |6|. 


IV. Results 

Now we exhibit some illustrative results demonstrating the 
developed algorithm’s performance. We investigate both the 
electric field £ z radiated by a vertical (i.e., ^-directed) Hertzian 
electric current dipole (VED), as well as the magnetic field H z 
radiated by a ^-directed Hertzian (equivalent) magnetic current 
dipole (VMD); both sources radiate at / = 2MHz. In both 
scenarios, the source resides at depth z’ = Om within a three- 
layer NBAM, occupying the region — 1 < 2 < 1 [m], of mate¬ 
rial properties e r = jl r = Diag[10,10,1/10], Diag[5, 5,1/5], 
and Diag[2, 2,1/2] within the regions — 1 < z < — 1/4 [m], 
— 1/4 < z < 1/4 [m], and 1/4 < z < 1 [m] (resp.); see 
Fig. [T] The top layer (z > lm) is vacuum (e r i = y r i = 1) 
while the bottom layer (z < — lm) is a perfect electric 
conductor (PEC); note that this layered-medium configuration 
was specifically chosen to facilitate comparison with closed- 
form solutions through invocation of T.O. and EM Image 
theory G3- Indeed the EM field solution within z > — lm, 
for our five-layered configuration involving a VED source, 
can be shown identical to the closed-form field result of two 
VED’s (located at depths 2 = —1.75m and 2 = —19.25m) 
of identical orientation to the original VED and radiating 
in homogeneous, unbounded vacuum. Note that within the 
NBAM, an added step to compute the closed-form result must 
be taken, appropriately mapping the observation points within 
the NBAM to vacuum observation points by viewing a d-meter 
thick NBAM layer e r = jl r = Diag[n, n, l/n\ as equivalent to 
a nd-meter thick vacuum layer. Similarly, the VMD problem 
can be shown identical to two VMD’s (located at depths 
z = —1.75m and z = —19.25m) radiating in homogeneous, 
unbounded vacuum; in this scenario however, image theory 
prescribes that the 2 = —1.75m VMD possess identical 
orientation to the original VMD, but that the z = —19.25m 
VMD possess opposite orientation Q 

Observing Figs. |2c||2d[ we note the relative errors in both 
the electric field (S e ) and magnetic field (Sh) are very low|^] 
approaching in most of the observation plane near the limits 
of floating point double precision-related numerical noise 
(approximately -150 in [dB] scale); for reference, Figs. |2a[ 
|2b| are the computed field distributions themselves from our 
algorithm. This is consistent with our having set an adaptive 
relative integration error tolerance of 1.2 x 10 -14 . We do 
observe however that the error noticeably increases (for fixed 
observer/source radial separation) as the observation angle 
tends closer to “horizon” (i.e., source depth z' and observer 

7 The amplitudes of the VED and VMD (i.e., lying within the central NBAM 
layer) must be scaled by a factor of 1/5 (relative to the vacuum sources) to 
facilitate field comparisons. Moreover the normal field components {£ z , H z }, 
within the NBAM layer with properties e r = p r = Diag[n, n, 1 /n\, are also 
scaled (artificially, for both visual display and error computation purposes) by 
1/n to account for their discontinuity across material interfaces. 

8 Let E c be the computed electric field, and let E v be the closed-form 
reference solution. Then S e = \E C — E V \/\E V \ (likewise for Sh)- 
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depth z coinciding). The error variation trend versus angle 
has been observed before urn even when the source resided 
in non-NBAM media, and hence the increased error versus 
observation angle is not likely due to instabilities in the 
presented NBAM-robust algorithm. We conjecture rather that 
the increasing error (versus observation angle) arises due to 
commensurately increasing numerical cancellatiorj^] that can 
only be partially offset by a (computer resource limited) finite 
extent of hp integration refinement performed using finite 
precision arithmetic. This numerical cancellation, we remark, 
is well known to be predominantly induced by integrand 
oscillation, which worsens as the observation angle tends to 
horizon (lj, p2| . One remedy is to use a constant-phase 
path 131, but a robust remedy for 2-D integrals (needed 
for generally anisotropic media) remains an open question. 
Moreover, this path would change as one varies the outer in¬ 
tegration variable. Finally we emphasize that given the design 
of our particular implementation, which always first computes 
the direct electric field and then (if need be) computes the 
magnetic field using ancillary relations jTTJ[Ch. 2], we have 
in fact tested the soundness of both Eqns. ( |III.l| )-( [Iir4| ) (VED 
scenario) and Eqns. ( |III.5| )-( |nL6| ) (VMD scenario). 

10m 



4m 


Fig. 1: Vertically-oriented Hertzian dipole current source within a 
three-layer NBAM. The purple (air) and blue (NBAM) regions form 
the plane on which the fields are observed in Fig. [2] The parameter 
n equals ten, five, and two within the regions —1 < z < —1/4 [m], 
— 1/4 < z < 1/4 [m], and 1/4 < z < 1 [m], respectively. 


V. Conclusion 

We addressed a fundamental origination of breakdown in the 
spectral-domain-based (PWE) evaluation of EM fields radiated 
by sources embedded within NBAM planar slabs, leading to 
a robust formulation that can accurately compute EM fields 
despite the modal degeneracy, induced by said NBAM, that 
would ordinarily lead to numerical instabilities and/or corrup¬ 
tion of the functional form of the plane wave eigenfunctions. 
Indeed this instability arises due to eigenvalues that, while non¬ 
defective, have an algebraic multiplicity equal to two rather 
than one. The remedy is to apply a proper (analytical) “pre¬ 
treatment” of the spectral-domain tensor operators prior to 

9 Namely, cancellation from radiation field contributions arising from nu¬ 
merical integration along contour sub-sections symmetrically located about 
the imaginary k x and k y axes. By contrast, our algorithm robustly ensures 
(irrespective of observation angle) that the evanescent field contribution intro¬ 
duces little numerical cancellation-induced error and rapid convergence 0. 


10Log 10 |E z | [dB] 10Log 10 |H z | [dB] 



-5 -4 -3 -2 -1 0 1 2 3 4 5 -5 -4 -3 -2 -1 0 1 2 3 4 5 


x-x' [m] x-x' [m] 

(c) (d) 

Fig. 2: (a) S z radiated by a VED. (b) Ji z radiated by a VMD. (c) 

Relative error: S z . (d) Relative error: Ji z . 

polarization amplitude extraction, resulting in robust analysis 
of EM fields in arbitrary anisotropic planar-layered media. 
Results validated the high accuracy of numerical computations 
based on this analytical pre-treatment. 
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